#include <stdlib.h>

#include "../Utils/Utils.h"
#include "../LocalStack/LocalStack.h"

#include "../SparceMatrix/sparse_matrix.h"
#include "../Containers/NumList.h"

#include "Vector.h"
#include "Matrix.h"

#include "ModularityMatrix.h"

elem* CreateModularityMatrix(const sparse_matrix_lst* adjMatrix, const NumList* group)
{
	elem* modulMatrix;
	
	START();
	
	Assert(adjMatrix != NULL);
	Assert(group != NULL);
	modulMatrix = (elem*)NEW(sizeof(elem) * group->size * group->size);
	modulMatrix = PartialModularityMatrix(adjMatrix, group, modulMatrix);
	modulMatrix = TopOfModularityMatrix(modulMatrix, group->size);

	RETURN(modulMatrix);
}

elem* PartialModularityMatrix(const sparse_matrix_lst* adjMatrix, const NumList* group, elem* modulMatrix)
{
	elem *degVector, invMatrixSum = 0;
	const NumList* groupIterI;
	const NumList* groupIterJ;
	cell_t* adjIter;
	unsigned modulIter = 0, modulInvIter = 0, rowCounter = 0;
	const unsigned modulSize = group->size;
	
	START();
	
	degVector = (elem*)NEW(sizeof(elem) * (adjMatrix->n));
	degVector = DegVector(adjMatrix, degVector);
	
	invMatrixSum = 1 / SumOfVector(degVector, adjMatrix->n);
	
	for (groupIterI = group; rowCounter < modulSize; rowCounter++, groupIterI = groupIterI->next)
	{
		modulIter += rowCounter;
		modulInvIter = modulIter - modulSize;
		
		adjIter = adjMatrix->rows[groupIterI->value];
		for (groupIterJ = groupIterI; groupIterJ->next != NULL; groupIterJ = groupIterJ->next)
		{
			while((adjIter != NULL) && (adjIter->colind < groupIterJ->value))
				adjIter = adjIter->row_next;
			
			modulMatrix[modulInvIter += modulSize] = (modulMatrix[modulIter++] = 
				((adjIter == NULL || adjIter->colind != groupIterJ->value)? 0 : 1) 
					- (degVector[groupIterJ->value] * degVector[groupIterI->value] * invMatrixSum));
		}
	}
	
	Delete(degVector);
	
	RETURN(modulMatrix);
}

elem* TopOfModularityMatrix(elem* matrix, unsigned size)
{	
	int iter = 0, row = 0;
	elem* vectorOfRowSums;
	
	vectorOfRowSums = (elem*)NEW(sizeof(elem) * size);
	vectorOfRowSums = VectorOfRowSums(matrix, vectorOfRowSums, size);
	
	for (iter = 0; row < size; iter += size + 1)
		matrix[iter] -= vectorOfRowSums[row++];
	
	Delete(vectorOfRowSums);
	
	RETURN(matrix);
}
